%% Ergasia1_splines_parembolh_12.03.2009
%% Onoma: Georgallides  Andreas    A.E.M : 723

%%%%%%%%%% Ergasia 1_a %%%%%%%%%%

clear
clc
disp('Kampylh Fysikhs Kyvikhs Paremvolhs');

%% INPUT
%% 1.
x= [1 2 3 4 5]   % meres
y= [7607.40 7544.60 7561.20 7526.40 7515.30];   % TSE300
z= [1848.20 1836.50 1860.20 1843.40 1805.20];   % NASDAC

n= length(y)-1;  %bathmos polyonymou = simeia - 1
h=1;   % 2-1, 3-2, 4-3, 5-4

%% 2.
%% Arxikopoihsh A Pinaka 
A = 4 * eye(n-1,n-1);
for  i = 1:n-1
     if(i == 1) 
         A(i,i+1)= 1;
     elseif(i == n-1)
         A(i,i-1) = 1;
     else
         A(i,i-1)= 1;
         A(i,i+1)= 1; 
     end
end
A
%% B dianysma gia y,z
for i=1:n-1
   by(i) =  (6/(h^2))* (y(i) - 2*y(i+1) + y(i+2)); % deikths TSE300
   bz(i) =  (6/(h^2))* (z(i) - 2*z(i+1) + z(i+2)); % deikths NASDAQ
end

%% 3.
%% Am=b  =>  m=A\b'
mtempy = A\by';  % TSE300
mtempz = A\bz';  % NASDAQ

% Ypologismos  m(j), j=1,2,...,n,n+1. gia TSE300, NASDAQ

for i=1:n+1
    if(i == 1)||(i == n+1)
        my(i) = 0;  % TSE300
        mz(i) = 0;  % NASDAQ
    else
        my(i)=mtempy(i-1);  % TSE300
        mz(i)=mtempz(i-1);  % NASDAQ
    end
end
my
mz


%% 4.
%% Ypologismos  paramertwn a(j), b(j), c(j), d(j) twn
%% polyonymwn kyvikhs paremvolhs y = Sj(x)

ay = y(1:n);
by = (y(2:n+1) - y(1:n))./h - h.*(my(2:n+1)+2*my(1:n))/6;
cy = 0.5*my(1:n);
dy = (my(2:n+1) - my(1:n))./(6*h);

az = z(1:n);
bz = (z(2:n+1) - z(1:n))./h - h.*(mz(2:n+1)+2*mz(1:n))/6;
cz = 0.5*mz(1:n);
dz = (mz(2:n+1) - mz(1:n))./(6*h);

%% 5.
%% Splines paremvolis sto diasthma [0,5] gia ta 2 synola

xInt = 1 : 0.5 : 5;

for j = 1 : length(xInt)
	if xInt(j) ~= x(n+1)
        iInt(j) = sum(x <= xInt(j));
	else
		iInt(j) = n;
	end
end
yInt = ay(iInt) + by(iInt).*(xInt-x(iInt)) + cy(iInt).*(xInt - x(iInt)).^2 + dy(iInt).*(xInt - x(iInt)).^3;
yEx = humps(xInt);
plot(xInt,yInt,'g',x,y,'b*');   axis([0,5,7500,7950]);

pause;
zInt = az(iInt) + bz(iInt).*(xInt-x(iInt)) + cz(iInt).*(xInt - x(iInt)).^2 + dz(iInt).*(xInt - x(iInt)).^3;
plot(xInt,zInt,'g',x,z,'b*');   axis([0,5,1800,1870]);

pause;


%%%%%%%%%% Ergasia 1_b %%%%%%%%%%
clear
clc
disp('Kyvikh spline paremvolhs me paravolikh simperifora sta akra');

%% INPUT
%% 1.
x= [1 2 3 4 5]   % meres
y= [7607.40 7544.60 7561.20 7526.40 7515.30];   % TSE300
z= [1848.20 1836.50 1860.20 1843.40 1805.20];   % NASDAC

n= length(y)-1;  %bathmos polyonymou = simeia - 1
h=1;   % 2-1, 3-2, 4-3, 5-4

%% 2.
%% Arxikopoihsh A Pinaka
A = 4 * eye(n-1,n-1);
for  i = 1:n-1
     if(i == 1) 
         A(i,i) = 5;
         A(i,i+1)= 1;
     elseif(i == n-1)
         A(i,i) = 5;
         A(i,i-1) = 1;
     else
         A(i,i-1)= 1;
         A(i,i+1)= 1; 
     end
end
A
%% B dianysma gia y,z
for i=1:n-1
   by(i) =  (6/(h^2))* (y(i) - 2*y(i+1) + y(i+2)); % deikths TSE300
   bz(i) =  (6/(h^2))* (z(i) - 2*z(i+1) + z(i+2)); % deikths NASDAQ
end

%% 3.
%% Am=b  =>  m=A\b'
mtempy = A\by';  % TSE300
mtempz = A\bz';  % NASDAQ

% Ypologismos  m(j), j=1,2,...,n,n+1. gia TSE300, NASDAQ
for i=2:n
        my(i)=mtempy(i-1);  % TSE300
        mz(i)=mtempz(i-1);  % NASDAQ
end
my(1) = my(2);  % TSE300
mz(1) = mz(2);  % NASDAQ

my(n+1) = my(n);  % TSE300
mz(n+1) = mz(n);  % NASDAQ
my
mz

%% 4. Ypologismos paramertwn  a(j), b(j), c(j), d(j) twn
%% polyonymwn kyvikhs paremvolhs y = Sj(x)

ay = y(1:n);
az = z(1:n);
   
by = (y(2:n+1) - y(1:n))./h - h.*(my(2:n+1)+2*my(1:n))/6;
bz = (z(2:n+1) - z(1:n))./h - h.*(mz(2:n+1)+2*mz(1:n))/6;

   
cy = 0.5*my(1:n);
cz = 0.5*mz(1:n);

dy = (my(2:n+1) - my(1:n))./(6*h);
dz = (mz(2:n+1) - mz(1:n))./(6*h);

%% 5. Splines paremvolis sto diasthma [0,5] gia ta 2 synola

xInt = 1 : 0.01 : 5;
for j = 1 : length(xInt)
	if xInt(j) ~= x(n+1)
        iInt(j) = sum(x <= xInt(j));
	else
		iInt(j) = n;
	end
end
yInt = ay(iInt) + by(iInt).*(xInt-x(iInt)) + cy(iInt).*(xInt - x(iInt)).^2 + dy(iInt).*(xInt - x(iInt)).^3;
yEx = humps(xInt);
plot(xInt,yInt,'g',x,y,'b*');   axis([0,5,7500,7950]);

pause;
zInt = az(iInt) + bz(iInt).*(xInt-x(iInt)) + cz(iInt).*(xInt - x(iInt)).^2 + dz(iInt).*(xInt - x(iInt)).^3;
plot(xInt,zInt,'g',x,z,'b*');   axis([0,5,1800,1870]);


pause;


%%%%%%%%%% Ergasia 1_c %%%%%%%%%%
clear
clc
disp('Kyvikh spline paremvolhs me eniaia kyvikh symperifora sta akra');

%% INPUT
%% 1.
x = [1 2 3 4 5];    % meres 
y = [7607.40 7544.60 7561.20 7526.40 7515.30];   % TSE300
z = [1848.20 1836.50 1860.20 1843.40 1805.20];   % NASDAQ

n = length(y)-1;  %bathmos polyonymou = simeia - 1
h = 1;  % 2-1, 3-2, 4-3, 5-4

%% 2.
%% Arxikopoihsh A Pinaka
A = 4 * eye(n-1,n-1);
for  i = 1:n-1
     if(i == 1) 
         A(i,i) = 6;
         A(i,i+1)= 1;
     elseif(i == n-1)
         A(i,i) = 6;
         A(i,i-1) = 1;
     else
         A(i,i-1)= 1;
         A(i,i+1)= 1; 
     end
end
A
%% B dianysma gia y,z
for i=1:n-1
   by(i) =  (6/(h^2))* (y(i) - 2*y(i+1) + y(i+2)); % deikths TSE300
   bz(i) =  (6/(h^2))* (z(i) - 2*z(i+1) + z(i+2)); % deikths NASDAQ
end

%% 3.
%% Am=b  =>  m=A\b'
mtempy = A\by';  %deikth TSE300
mtempz = A\bz';  %deikth NASDAQ

% Ypologismos  m(j), j=1,2,...,n,n+1. gia TSE300, NASDAQ
for i=2:n
        my(i)=mtempy(i-1);  % TSE300
        mz(i)=mtempz(i-1);  % NASDAQ
end
my(1) = (2*my(2)) - my(3);     % TSE300
my(n+1) = (2*my(n)) - my(n-1); % TSE300
mz(1) = (2*mz(2)) - mz(3);     % NASDAQ
mz(n+1) = (2*mz(n)) - mz(n-1); % NASDAQ
my
mz

%% 4. Ypologismos paramertwn  a(j), b(j), c(j), d(j) twn
%% polyonymwn kyvikhs paremvolhs y = Sj(x)
ay = y(1:n);
by = (y(2:n+1) - y(1:n))./h - h.*(my(2:n+1)+2*my(1:n))/6;
cy = 0.5*my(1:n);
dy = (my(2:n+1) - my(1:n))./(6*h);

az = z(1:n);
bz = (z(2:n+1) - z(1:n))./h - h.*(mz(2:n+1)+2*mz(1:n))/6;
cz = 0.5*mz(1:n);
dz = (mz(2:n+1) - mz(1:n))./(6*h);

%% 5. Splines paremvolis sto diasthma [0,5] gia ta 2 synola
xInt = 1 : 0.01 : 5;
for j = 1 : length(xInt)
	if xInt(j) ~= x(n+1)
        iInt(j) = sum(x <= xInt(j));
	else
		iInt(j) = n;
	end
end
yInt = ay(iInt) + by(iInt).*(xInt-x(iInt)) + cy(iInt).*(xInt - x(iInt)).^2 + dy(iInt).*(xInt - x(iInt)).^3;
yEx = humps(xInt);
plot(xInt,yInt,'g',x,y,'b*');   axis([0,5,7500,7950]);

pause;
zInt = az(iInt) + bz(iInt).*(xInt-x(iInt)) + cz(iInt).*(xInt - x(iInt)).^2 + dz(iInt).*(xInt - x(iInt)).^3;
plot(xInt,zInt,'g',x,z,'b*');   axis([0,5,1800,1870]);
